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We present a numerical code for calculating the self force on a scalar charge moving in a bound 
(eccentric) geodesic in the equatorial plane of a Kerr black hole. We work in the frequency do- 
main and make use of the method of extended homogeneous solutions [Phys. Rev. D 78, 084021 
(2008)], in conjunction with mode-sum regularization. Our work is part of a program to develop a 
computational architecture for fast and efficient self-force calculations, alternative to time-domain 
methods. We find that our frequency-domain method outperforms existing time-domain schemes 
for small eccentricities, and, remarkably, remains competitive up to eccentricities as high as ~ 0.7. 
As an application of our code we (i) compute the conservative scalar-field self-force correction to the 
innermost stable circular equatorial orbit, as a function of the Kerr spin parameter; and (ii) calculate 
the variation in the rest mass of the scalar particle along the orbit, caused by the component of the 
self force tangent to the four-velocity. 

I. INTRODUCTION 

The last few years have seen a breakthrough in the program to compute the gravitational self-force (SF) on massive 
particles in orbit around black holes (see, e.g., jl| and references therein). Such systems are employed as models for 
strongly gravitating astrophysical binaries of extreme mass ratios (compact stellar objects captured by massive black 
holes in galactic nuclei), which are key sources of gravitational waves for planned detectors. The long-term ambition 
of the SF program is to model the phase evolution of the inspiral system, and the emitted gravitational waveforms, for 
generic inspiral orbits about rotating (Kerr type) black holes. However, to date, all computations of the gravitational 
SF have assumed a non-rotating (Schwarzschild type) central hole. The generalization to Kerr spacetime presents a 
major technical challenge, which the community must now come to address. 

In a recent paper [2j (hereafter "Paper I" ) we reported a first computation of the SF for an orbit around a Kerr black 
hole. In this computation we resorted, as often in the SF program, to the simple framework of a scalar-charge toy 
model, and considered the scalar-field SF (SSF) acting on such a particle as it moves along a (fixed) circular equatorial 
geodesic of the Kerr geometry. The scalar-charge model provides a convenient environment for development, while 
already capturing much of the technical complexity of the full gravitational problem. 

In Paper I we chose to take a frequency domain (FD) approach to the SF problem: The (scalar) field equations 
are decomposed into Fourier-harmonic modes on the Kerr background, and the resulting fully-separated ordinary 
differential equations (ODEs) are integrated numerically with suitable boundary conditions, in what is a routine 
procedure in black-hole perturbation studies. The Fourier-harmonic modes are then post-processed and used as input 
for the standard mode-sum regularization formula Q that yields the SSF. The FD approach has the obvious advantage 
that its numerical component only involves the solution of ODEs. As we have seen in Paper I, this makes it extremely 
efficient computationally in comparison with existing time-domain (TD) algorithms, which involve the numerical 
evolution of partial differential equations (PDEs). In the circular-orbit case, typical run-times of a current-day TD 
code (to compute the SF at a single radius) are of order hours, whereas an FD code arrives at the answer in seconds. 
The performance of the FD method is expected to deteriorate with increasing orbital eccentricity, as the Fourier 
spectrum of the wave equation broadens [H, but, as we shall see in the current work, the method remains an efficient 
alternative to TD even for eccentricities as large as ~ 0.7. 

Two main potential drawbacks of the FD approach, in view of the ultimate goals of the SF program, are the 
following. First, it is not obvious how an FD scheme might be implemented in an efficient way in a self-consistent 
evolution scheme that solves the field equations simultaneously with the equation of motion in order to track the self- 
forced evolution of the orbit. While TD algorithms can, in principle, feed the SF information back into the equation 
of motion "in real time" , it is yet to be seen how FD schemes could accommodate a slowly evolving spectrum. One 
could possibly use FD data within the "osculating geodesies" approach to the orbital evolution (see [1, @), which 
requires as input the value of the SF along each member of a sufficiently dense sequence of fixed geodesies. At any 
rate, FD codes can provide an extremely useful machinery for testing the results of TD algorithms. They can be used 
to produce a large amount of data quickly and accurately, and the core methods they rely on are entirely independent 
from those of the TD codes. 

The second potential drawback becomes apparent when one moves on to consider the gravitational SF problem on 
the Kerr background. To the best of our knowledge, there is no known way to fully decompose the Lorenz-gauge 
perturbation equations in Kerr into Fourier-harmonic modes. It might be possible to proceed by using the tensorial 
spherical-harmonic decomposition as in the Schwarzschild problem [J 043 ano - then properly account for the coupling 
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between different harmonics that would occur in the Kerr case. An alternative would be to follow the radiation-gauge 
strategy advocated by Friedman et al. [H|[Hj]> which is naturally designed for an FD treatment. 

In the current work we remain focused on the scalar field problem, and present an extension of the analysis of Paper 
I from circular (equatorial) orbits to eccentric (equatorial) orbits. This generalization is highly non-trivial in an FD 
treatment. The move from a single fundamental frequency (in the circular orbit case) to a bi-periodic motion brings 
with it several new technical complications (see below), which must be addressed. We develop here computational 
tools to deal with these complications; these tools should be transferable to the gravitational problem. This work 
represents a first implementation of the standard mode-sum regularization technique [12L Il3j for non-circular orbits in 
Kerr. As such, it contains, as a secondary result, a first numerical confirmation of the regularization parameter values 
for such orbits (these analytic parameters were derived long ago Q but were never used in an actual computation 
thus far). 

Eccentric-orbit calculations in the FD, in both Schwarzschild and Kerr spacetimes, are hampered by the poor 
convergence of the Fourier sum near the particle singularity. The mode-sum formula requires as input the spherical- 
harmonic Z-modes of the retarded field (and their derivatives), evaluated at the particle and summed over all uj- 
frequencies. A naive attempt to reconstruct the Z-mode derivatives — which are generally discontinuous at the particle — 
as a sum of smooth Fourier harmonics leads to Gibbs-type oscillations, which severely impair the convergence of the 
Fourier sum in the vicinity of the particle (and, in fact, yields the wrong value at the particle itself). This problem has 
been identified and analyzed in Ref. [Tij , where a simple solution that entirely circumvents the problem was proposed. 
Using the proposed technique, dubbed the method of extended homogeneous solutions (EHS), the physical (retarded) 
field and its derivatives are reconstructed as an exponentially convergent sum over certain homogeneous solutions of 
the field equation. This technique has already been applied successfully in several numerical studies of perturbations 
from point particles 0, dUdEI, but here we implement it for the first time in a full computation of the SF. For that 
reason, we will take some time, in Sec. IVI C[ to assess the computational performance of EHS. 

A second challenge in extending the analysis of paper I to eccentric orbits has to do with a certain feature of 
the mode-sum scheme already mentioned above: the scheme requires as input the /-modes of the field (and its 
derivatives) in a sp/ierzcaZ-harmonic decomposition. Unfortunately, the field equations in Kerr spacetime do not 
separate into spherical harmonics, instead separating into spheroidal-harmonic modes (in the FD). This means that 
the numerically-computed spheroidal modes must be post-processed in a procedure that projects them onto a basis 
of spherical harmonics. In paper I we demonstrated the manageability of such a procedure in the circular-orbit case, 
where the spheroidicities encountered (proportional to the square of the modal frequencies) are relatively small, leading 
to a relatively weak coupling between the spheroidal and spherical modes. A potential concern is that the coupling 
might become less tractable with increasing eccentricity, since high harmonics of the orbital frequencies become more 
important, and these have higher spheroidicities leading to stronger coupling. We show here, however, that even at 
moderately large eccentricities the coupling remains weak enough to remain manageable in practice. 

Calculations of the SF for eccentric orbits give access to some interesting information about the post-geodesic 
dynamics in the binary system. In the gravitational problem, SF results have recently been used to quantify the 
precession effects of the conservative piece of the SF, and derive the resulting shift in the location and frequency of 
the innermost stable circular orbit (ISCO) I n the SSF problem, the ISCO shift in Schwarzschild spacetime 

was derived by Diaz-Rivera et al. Here we shall use our eccentric-orbit code to explore two effects of the SSF. 
First, we shall extend the work of Ref. to the Kerr case, and compute the conservative shift in the location and 
frequency of the ISCO (for equatorial orbits) as a function of the Kerr spin parameter a. Second, we will compute the 
variation in the particle's rest mass as it moves along the eccentric orbit. This variation is caused by the component 
of the SSF tangent to the particle's four velocity. We will verify that this effect does not lead to a net change in the 
rest mass over a full orbital revolution, as expected for the configuration at hand. 

The remainder of this paper is structured as follows. In section [TT] we review the relevant features of eccentric 
equatorial geodesies in Kerr geometry, and describe the formalism governing scalar-field perturbations. In section 
IIIII we discuss the application of the mode-sum scheme to orbits in Kerr spacetime. In section IIVI we examine the 
difficulties encountered with the naive FD approach to SF calculations for eccentric orbits and review the resolution 
to these problems using the method of EHS. In section [V] we provide details of our numerical scheme and in section 
IVII we present a sample of our results and discuss the computational performance of the method. Sections I VIII and 
IVIIII contain analyses of the conservative ISCO shift and of the rest- mass variation, respectively. Lastly in section ITXl 
we summarize our results and consider future work. Throughout this paper we use Boyer-Lindquist (BL) coordinates 
(t, r, ip) with metric signature ( — h H — h) and geometrized units such that the gravitational constant G and the speed 
of light c are equal to unity. 
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II. PRELIMINARIES 
A. Orbital setup and equation of motion 

Consider a pointlike particle of mass /z carrying a small quantity of scalar charge q and moving on a bound orbit 
about a Kerr black hole of mass M 3> /i and spin aM. We assume that the correction to the background Kerr 
spacetime from the stress-energy of the particle's scalar field can be neglected. We denote the particle's worldline 
(in Boyer-Lindquist coordinates) by x£(r) and its four- velocity by u a (r) — dx a (r)/dT^ where r is the proper time. 
In this work we ignore the gravitational SF and seek to calculate only the SSF [of 0(q 2 )}. The orbital dynamics is 
described by the (self-)forced equation of motion [2(| 

uPVpQnf) = qV a $ R = F» 1£ , (1) 

where is Detweiler- Whiting's smooth, regularized field [2l| . and the covariant derivative is taken with respect to 
the background Kerr metric. In this work we calculate the SSF F s " lf (oc q 2 ) that would be felt by a particle moving 
on a fixed geodesic orbit, envisaging that this information could later be used to compute the true inspiralling orbit 
using a scheme similar to that described in In this work we only consider the case of bound equatorial (8 p = tt/2) 
orbits. 

As discussed in Paper I, once initial conditions are specified, a geodesic orbit about a Kerr black hole is uniquely 
parametrized by three constants of motion: its specific energy £ = — Ut, angular momentum C = u v and Carter 
constant Q. In the case of equatorial orbits the Carter constant vanishes and thus the pair (£,C) suffices to specify 
the orbit. It turns out to be convenient to use a different pair of parameters consisting of the semi-latus rectum p and 
eccentricity e, representing strong-field generalizations of their Keplerian counterparts. For eccentric orbits we denote 
the BL radius at the point of closest approach (periastron) and the BL radius when the two bodies are furthest apart 
(apastron) by r m - ln and r max , respectively. Then p and e are defined by 

2^max^min ^max ^*min /~\ 

e = z ~z — • ( 2 ) 



+ r n 



The relation between (p, e) and the specific energy and angular momentum is found to be [2^ 
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C = x + a£ , (4) 

where x = x{a 1 p,e) is a rather complicated function given in appendix lAl All bound (stable) equatorial geodesies 
havep 2 > a; 2 (l + e)(3 — e); the (a dependent) curve p 2 = x 2 (l + e)(3 — e) defines a separatrix in the e-p plane (see, e.g., 
Fig. 2 of [13]), with all orbits below the separatrix being unstable. For each given a, the intersection of the separatrix 
and the axis e = defines the ISCO. 

Using the p, e paramctrization, the particle's orbital radius in BL coordinates can be expressed as 

p 

r = r p(x) = -r l > (s) 

1 + e cos x 

where \ is a monotonically-increasing parameter along the particle's worldline. The azimuthal angle ip p and coordinate 
time t p can be computed as functions of x along the particle's worldline using the expressions given in appendix [X] 
We take <p p to be monotonically increasing for all orbits (or, equivalently, £ > 0) and distinguish between prograde 
and retrograde orbits by the sign of a (a > for prograde, a < for retrograde). Without loss of generality we 
let t p (x = 0) be the time of a periastron passage and denote the radial period (the t-time taken for the particle to 
progress from a periastron passage to a subsequent one) by T r = t p (\ — 2tt) = 2t p {\ = vr). Similarly we denote the 
change in cp p over a time period T r by A.ip p = (p p (x = 27r) = 2(p p (x = 7r). Then the two frequencies associated with 
the radial and azimuthal motion are given by 

<V = £, (6) 



Note that in Eq. (Jl} we have kept the mass term \x inside the derivative operator. Expanding the derivative one finds 
a term orthogonal to the particle's four- velocity, which is responsible for the self-acceleration, and a term tangential to 
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the four-velocity, which in general gives rise to a dynamically varying mass. The orthogonal and tangential components 
of Eq. (P) are given respectively by 

dn a 

M^r = ( S P + u"«f>)FL = F±W > (7) 
^ = -u a Ff { . (8) 

Combining Eqs. ((TJ and ([8} allows us to write the mass change explicitly as a function of r, 

M (r) - Mo - q<f> R (r) , (9) 

where /x is a constant of integration (sometimes called the bare mass). In the case of circular, equatorial orbits 
discussed in Paper I there was no change in the particle's rest mass due to the stationarity of the setup. For eccentric 
orbits this is no longer the case. However, since $> r (t) comes back to itself after a radial period T r , we expect no net 
mass change over that period. 



B. Field equation and multipole decomposition 

We assume that the field $ associated with the particle's scalar charge q obeys the minimally coupled Klein-Gordon 
equation 

V Q V Q $ = -4ttT , (10) 
where T denotes the particle's scalar charge density. We model the latter as 

T = q [ *V - x^l-gix)- 1 ' 2 ]^ = J-S(r - r p )6(v> - <p p )6(6 - tt/2) , (11) 

where g = —p 4 sin 2 9 (with p 2 = r 2 + a 2 cos 2 9) is the determinant of the background Kerr metric, and in the second 
equality we have specialized to equatorial orbits with 9 p = tt/2. The t component of the four-velocity, it', is calculated 
as u l — g tip £ — g tt £ 1 where g a/3 are the contravariant components of the Kerr metric tensor in BL coordinates, 
evaluated at the particle. 

The scalar wave equation (|10p in Kerr geometry can be completely separated into spheroidal harmonic and frequency 
modes in the form [23|, HiJ 

« oo I 

* = J E E R imJr)S im (9; a 2 )e^e-^ . (12) 
f=0 m=-I 

Here S~ lm (9\a 2 ) are spheroidal Legendre functions [see Eq. (|2U|) below] with spheroidicity 

(j 2 = -a 2 uj 2 . (13) 

We reserve the term spheroidal harmonic for the product S^ m {9; a 2 )e tmv . Notice that we label spheroidal-harmonic 
modes by Im as we will later introduce sp/ierzcaZ-harmonics modes which we will label by Im. The spheroidal harmonics 
are orthonormal with normalization given by 

S tm (9; o- 2 y mv S Pm , (9- o- 2 )e~ im '* dQ = 8 n ,5 mm , , (14) 

where d£l = sin 9d9dip, and 5 ni „ 2 is the standard Kronecker delta. 
The source T in Eq. (1111) has a discrete spectrum given by 

u = uj mn = mO v + n£l r , (15) 

for integer m and n. It therefore admits a discrete Fourier decomposition, of the form 

oo / oo 

A-EE E T u(«^"^'. ( 16 ) 

f=0 m=-i n =-°° 
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where Sj mn (9) = St (6; — a 2 ^ 2 ^), and the factor p 2 is introduced for later convenience. Using the orthonormality 
property (|14[) and taking the inverse Fourier transform of Eq. (|16p we find 

T tmn (r) = qS irnn ^/2)T^ r^y^lr - r p {t)]e^~ m ^t)\ dt , (17) 

Jo 

where = ip p (x(t)), and x(<) is obtained formally by inverting t p {x) (this inverse exists since t is a monotonically 

increasing function of x). Noting that TV (r) only has support on r min < r < r max and changing the integration 
variable from t to r p [taking t p {\ = 0) = <p p (x = 0) = without loss of generality ], we finally obtain 

T tmn (r) = 2g ^p ( ( ^ 2) cos[uW P (r) - m^ p (r)} x 6(r - r min ) x 9(r max - r) , (18) 

where is the standard Heaviside step function, and t p (r) and tp P {r) are obtained by formally inverting r p (x) hi the 
range < x < 7r. 

The decompositions (fT6]l and (the discrete version of) (|T2l) separate the r and dependence of the field equation 
(1101) . with resulting (mn-mode equations given by 

^d~r ( A ^ HL j + [fl2m2 _ 4M ~™« + + a2 f"ln ~ aW mn A - A fm A)] fl fm „ = -4,AT im Jr) , (19) 
1 d ( . „dSt \ /, , , 9 „ TO 2 



sin 6*961 V d6> 



sin ) + ( x tm + a 2 u 2 mn cos 2 6--^)S tmn = 0, (20) 



r" — 2Mr + a 2 and R~ lmn = R{ mul ■ The angular equation (|2"0)) takes the form of the spheroidal Legendre 
equation with spheroidicity a 2 = —a 2 ^^. Its eigenfunctions are the (frequency dependent) spheroidal Legendre 
functions Si mn (9;u 2 ), with corresponding eigenvalues A| . In general there are no closed-form expressions for Sj mn 
or A ; ~ m but they can be calculated using the spherical-harmonic decomposition method described in appendix A of 
Paper I. For a 2 = the spheroidal harmonics Si mn e %mip become the standard spherical harmonics Yt , and their 
eigenvalues reduce to Xt = Z(Z + 1). 

To simplify the construction of boundary conditions (below), and assist our numerical scheme, it is convenient to 
transform to a new radial variable, 

= rR imn (r) , (21) 
and introduce the tortoise radial coordinate r* defined through 

dr A < 22 > 
With the above definition the tortoise coordinate is given explicitly in terms of r as 

- = + M '"<^' 2 > + '» ■ < 23 > 



where we have specified the constant of integration, and where r± = M ± V M 2 — a 2 are the outer and inner roots 
respectively of the equation A = 0. There are several standard alternative choices that can be made for the tortoise 
coordinate. Our choice is motivated by the observation of Bardeen et al. [25[ that it leads to a simpler radial potential 
than other common choices [24j |. This, in particular, simplifies the construction of the numerical boundary conditions 
for the resulting radial equation. In terms of ipt m ( r ) an d t* , the radial equation (|19l) takes the simpler form 

+ V tmn (r)i> tmn = -^T tmn (r) = Z imn (r) , (24) 
where Tt is given in Eq. (|18p above, and Vt is an effective (w-dependent) radial potential given by 



V t m J r ) 



(r 2 + a 2 )uj mn — ami 2 A [ 2 2 2(Mr — a 2 



A 

74 



Im 



2amuj mn + a ui mn + 



(25) 



There is no known closed-form analytic solution to the radial equation (|24[) for general Imn and so it has to be solved 
numerically. Our technique for doing this is detailed in section [V] 
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C. Boundary conditions 

The solution of the radial equation (|24l) is uniquely determined outside the black hole once suitable boundary 
conditions are specified on the horizon (r* — > — oo) and at spatial infinity (r* — > oo). In order to select the retarded 
solution to the field equation (ITU1) we require that radiation be 'outgoing' at infinity and 'ingoing' at the horizon. 
Making this statement precise is straightforward at spatial infinity but there is some subtlety at the event horizon 
following from the dragging of inertial frames in Kerr spacetime. We go though the calculation for both boundaries in 
section II. C. of Paper I and as the boundary conditions have no dependence on the type of orbit being examined all 
that is said there carries through here, with the only change being that uj is now bi-harmonic in to and n. Here it will 
suffice to simply state the results presented in Paper I. The asymptotic boundary conditions for the radial equation 
at spacial infinity and the event horizon are given, respectively, by 

Vw(r* -> oo) ~ e +iu ™ r * , (26) 
Vw(r* -> -oo) ~ e -*iw. , (27) 



where 



2Mr + uj mn - am 

7mn = 2 ■ ^o) 



These boundary conditions specify a ratio between the value of the radial function and its derivative. Later on, in 
section IV A[ we will set an (arbitrary) amplitude scaling for the homogeneous solutions. The unique inhomogeneous 
solutions will then be constructed via the method of EHS, as described in section ITV Al 



III. SELF-FORCE VIA MODE-SUM REGULARIZATION 



In the mode-sum scheme the full field is decomposed into spherical-harmonic modes, even in Kerr spacetime. We 
define the full force at an arbitrary field point x = x^ in the neighbourhood of our scalar charge q as the field 

Ff\x) ee qV a $(x) = Fi ml) \x) > (29) 
l 

where $(a;) is the physical (retarded) scalar field, and i^ fu11 ^ denotes the total contribution to gV a $ from its spherical 
harmonic /-mode (summed over to). Each F^^' 1 is finite at the particle's location, although in general the sided 
limits r — ¥ r^ yield^two different values, denoted F^±^ 1 respectively. The mode-sum formula for the SSF acting on 



the particle reads [1] 



1=0 



(f'± u11) -A a± (l + 1/2)- B a ) , (30) 



were the coefficients A a ± and B a are /-independent "regularization parameters" , the values of which are known 
analytically for generic orbits about a Kerr black hole [3[ (see Ref. [lj] for a full derivation of the regularization 

parameters in the Kerr case). We note the difference F^J .(x p ) — A a ±(l + 1/2) no longer exhibits the ± ambiguity. 

In Kerr, as we have seen, the field naturally decomposes into spheroidal harmonic modes. The mode-sum scheme, 
on the other hand, requires spherical harmonic modes as input even in Kerr spacetime. Hence in order to regularize 
using the standard mode-sum approach we first need to project the spheroidal-harmonic contributions onto a basis of 
spherical harmonics. We do this by expanding 

oo 

S tm {6; a 2 K mv = E &L(O*W0> f) , (31) 

1=0 

where the cr-dependent coefficients b\ m are determined from a recursion relation found by substituting the series expan- 
sion into the angular differential equation ([20]) — a procedure which follows [26[ and is described in detail in Appendix 
A of Paper I. Given the coupling coefficients b l lm , we can write the spherical-harmonic ("Z-mode") contribution to the 
full force as 



E <f>im(t,r)Yi m (9,(p)/r 



.m=—l 



(32) 
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where 

oo oo 

<t>lm(t, = E E bim^iJ^e-^ . (33) 
«=°i=|m| 

In the last equation the convergence of the sum over n is impaired near the particle due to the Gibbs phenomenon; 
we will discuss this problem and its resolution in section ITVl below. 

Formally, when constructing <j)i m via Eq. (|33|) one has to sum over an infinite number of spheroidal I modes. In 
practice this is not necessary as the coupling between the spheroidal and spherical harmonic modes is relatively 
narrow-band for the spheroidicities encountered in our calculation. In Paper I we numerically demonstrated that the 
contribution from a given spheroidal Im mode to the spherical harmonic Im modes of the field is strongly peaked around 
I = I with exponential wings. The coupling strengthens as the magnitude of the spheroidicity (= a 2 ui 2 ) increases. In 
the case of circular orbits we found that for the largest spheroidicity encountered in our work (for a black hole spin of 
a = 0.998A/ and orbital radius of ro = 2M) the coupling was still weak enough to be perfectly manageable in practice. 
For example, computing the first 50 spherical-harmonic mode contributions required only ^56 spheroidal-harmonic 
modes (i.e. even in this extreme case the bandwidth of the coupling is still only ±6 modes). In the case of eccentric 
equatorial orbits the situation is slightly worse as the spectrum is now bi-periodic (uj — m£l v + n£l r ), and the relevant 
values of uj, and consequently a 2 , can be much larger than in the circular orbit case. Nonetheless we find that the 
coupling is still weak enough to allow for reasonably efficient SSF computations (see section IVI CI for details of the 
computational performance of the method). As an example, for an orbit with parameters (a,p, e) = (0.9M, 10M, 0.5) 
(giving the highest spheroidicity considered in this work) we find, for / = 14, an effective bandwidth of ^12 modes. 



A. Conservative and dissipative pieces of the SF 

The SF can be split into conservative and dissipative components [IH [27j • This split is useful for analysing the 
different physical effects of the SF. Splitting the SF into its conservative and dissipative components is also practically 
beneficial as the two pieces admit Z-mode sums with different convergence properties, which are better dealt with 
separately. Let us write 

pself _ pret peons _i_ pdiss (^A^ 

r a — r a ol ' ol 1 

where 

F™™ = \ {FT + ^ Q adv ) , Ff™ = \ [FT ~ F^) , (35) 
in which the advanced SF is obtained from Eq. (|3Q|) by replacing Fa full ^($ r et) with F^ fu11 ^ ($ a dv), i-e. 

oo 

F adv = ^ [F ^»($ adv ) - A a± (l + 1/2) - B a ] , (36) 
z=o 



where ^adv is the advanced solution to Eq. (|TC 

Substituting into formulae (|35[) from Eqs. (|30[) and (|36p . we obtain mode-sum regularization formulae for the 
conservative and dissipative components [13j : 



E [ x o' 

(=0 



F ?un(con S) _ A ± {l + 1/2) _ B ^ (37) 

F d iss = J2F^l {diss) , (38) 



^,full(diss 

where 

1 rP fu11 {(b \ 4- F fu11 (cb M pfull(diss) _ 1 



F f u n ( co„ s) ^ i ($rct) + piun ($adv)]j F fun ( d iss ) ^ _ ^($ adv )] . ( 39 ) 



Note that the dissipative piece of the SF does not require regularization, whilst the conservative piece does. It 
turns out that the l-mode sum for the dissipative piece converges exponentially, while that for the conservative piece 
converges only as oc l^ 1 (the l-mode terms fall off as cx l/l 2 at large I) (l3l |. 
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The advanced field $ a dv can be constructed by reversing the boundary conditions (l26l) and ([27)) so that radiation is 
ingoing at infinity and outgoing at the horizon. This, however, doubles the computational cost. Fortunately there is a 
simple way around this, which takes adva ntag e of the particular symmetries of geodesies in Kerr (noted in a different 
context by Mino in (28[): As explained in [l3j it is possible to re-express Eq. (|55|) in terms of the retarded force F™* 
alone, by considering the value of F^ et at two "opposite" points along the equatorial orbit, i.e., ones with the same 
value of r p but opposite radial velocities. It is convenient to take r = to correspond to a periapsis passage, and 
express the SSF as a function of r along the orbit. Then one can show [l3| 



^ on V) - ^K C 'M + ewCH], ^ Q diss (r) = -K et W - H«)K ct (-T)} , (40) 



\[K ct (r) + £(a) Jf(-r)], Ft s (r) = \ 

where £( Q ) = (—1,1,1,-1) in Boyer-Lindquist coordinates. Analogous relations apply for i^"^ co and 8 
and we shall use them in our analysis to construct the necessary input for the mode sums (|37l) and (1381) . 



IV. REVIEW OF THE HIGH-FREQUENCY PROBLEM AND ITS RESOLUTION 

This problem was first discussed in Ref. [14| . We review it here for completeness, and where appropriate generalize 
the discussion to Kerr spacetime. 

Let ijjf (r) and i/j7 (r) De t w0 homogeneous solutions to the radial equation (|2"4"|) . satisfying the boundary condi- 
tions presented in Eqs. ([25]) and ([27]). respectively. These two homogeneous solutions form a basis that can be used 
to construct the physical inhomogeneous solution using the standard variation of parameters method: 

r r ib7 (r')Z; (r')r' 2 f»« i/>+ (r')Z; (r')r' 2 
jh. (r) - xh + (r) / lmn dr' + %b7 (r) / lmn dr' 

- CM • (41) 
Here we have changed the integration variable from r* to r using Eq. (1221) . and 

W = tp7 (dipt /dr*) - ipf (dip7 /dr*) = const (42) 

lmn lmn lmn lmn 

is the Wronskian. In the regions r < r m i n and r > r max this formula reduces to the homogeneous solutions 

f \ \ C 7 ^7 = $7 ( r )' r - r min , 

7 /,, ( r \ — J lmn lmn lmn ( 

Vw(r)-< + + - + r > w , ^ 

V lmn lmn lmn 

where the coefficients C~ and are given by 

lmn lmn 

f r ~ tbf <r)Z; (r)r 2 
c ± =w -i ZJmnl 1 mnK ' dr . (44) 

lmn J r ^ A(r) V ; 

The source term Zf (r), defined in Eq. (|24l) . has singularities at r — r max , r m j n and to avoid them it is convenient to 
change the integration variable from r to t p (r). Taking t p (r m i n ) = 0, the scaling coefficients Cf 1 are then given by 

xk 87r lS imn (7r/2) f T ^ 2 ilif m Jr p (t))cos(L} mn t-m<p p (t)) 



l ™~ T r W J r p {t)u\r p {t)) ' W 

where now the integrand is free from singularities. 

We now need to construct the time-domain field <fri m (t,r) and its t and r derivatives. Recalling Eq. (1331) . this field 
is constructed via 

<Pim(t,r) = ^ (f>lmn(t,r) , (46) 
n 

with 

oo 

<t>lmn(t,r) = £ 6 l^M e "" m "*' (47) 

[=0 
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where the b l lm 's are the spheroidal-spherical harmonic coupling coefficients from Eq. (|33~|) - For any given ro between 
r m j n and r max , the field <f>i m (t, ro) is not a smooth function of time (its derivative suffers a jump discontinuity when the 
particle crosses ro). Standard Fourier theory tells us that the sum over n in Eq. (l46l) will suffer Gibbs-type oscillations 
near the particle, and will converge very slowly as a result. Even more troubling, the derivatives of the field, <j>i m a , 
which are needed as input for the mode-sum scheme, may fail to converge to the correct value at the discontinuity. 
One might attempt to extract the correct value of the derivatives through extrapolation, but such a procedure would 
be computationally inefficient because of the poor convergence of the n-mode sum in the vicinity of the particle. 

It should be clarified that the above problem does not occur in the case of circular orbits, even if they are inclined, 
as for these orbits the field at any fixed radius is a smooth function of time. 



A. Method of extended homogeneous solutions 

The above problem might have made SF calculations in the frequency domain rather unattractive. Fortunately, Ref. 



14j recently proposed a simple technique for overcoming this problem, dubbed the method of extended homogeneous 
solutions (BBS). Ref. Q outlined the details of the method and provided a numerical example using EHS to calculate 
the monopole contribution to the scalar field for a particle in an eccentric orbit about a Schwarzschild black hole. 
Here we will overview the method and provide the formula required to extend their calculation to eccentric equatorial 
orbits about a Kerr black hole. 

The first step in the method (as its name suggests) is to consider an extension of the homogeneous solutions ib^ 

Imui 

to the entire domain, defined through 

4>f (r) = Cf 1>? (r), r>2M, (48) 

Iran ' Iran Iran ' v 1 

where the coefficients and C~ are those given in Eq. (44). Note that ibf and tl>7 coincide with in the 

Iran Iran Iran Iran Iran 

respective domains r > r max and r < r m i n , but neither is a solution of the inhomogeneous n-mode equation (|24[) in the 
sourced domain r m - m < r < r max . One then defines the two spherical-harmonic time-domain extended homogeneous 
solutions <t>f m and (j)~ n by 

tL(t,r)^J2tLn(t>r), (49) 

n 

where 

oo 

f=o 

Ref. [IH showed that the n-mode sum in Eq. (|4"!)|) converges exponentially-fast in |n|, and that it does so uniformly in 
t and r throughout r > 2M (note that (j>f m are smooth functions of t for any 7-). Furthermore, Ref. [l4| argued that 
the extended homogeneous solutions can be used to construct the true time-domain function <pi m {t,r) on either side 
of the particle's trajectory: 

The SSF is then obtained via Eqs. and (|30p as usual. 

The advantages of this method are clear. First, the problem of Gibbs ringing encountered when constructing 
the derivative of the field is no longer present; the derivative of the field converges with n to the correct value 
everywhere, including at the location of the particle. Second, the convergence of both the field and its derivatives is 
exponential, making an EHS calculation extremely computationally efficient. Third, the standard Wronskian-based 
method requires the evaluation of two integrals [see Eq. (|41[) ] for each orbital point separately (and for each Imn) . In 
EHS, on the other hand, the two integrals in Eq. (14"4"|) suffice to determine the Imn-mode along the entire orbit. 

Since the introduction of the method of EHS it has been applied to a range of problems, including the calculation 
of the monopole and dipole contributions to the gravitational SF in the Lorenz gauge [1], and the calculation of the 
full gravitational perturbation in the Regge- Wheeler gauge . Our work represents a first implementation of EHS 
for a full calculation of the SF. 
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V. NUMERICAL IMPLEMENTATION 

The radial equation (|24l) has no known closed-form analytic solution for general Imn and as such it needs to be 
solved numerically. For a given point in the (a,p, e) parameter space, a calculation of the SSF will entail solving for 
many imn modes (see below for details). To assist the calculation we note that ipf remains unchanged under the 
combined operations of (m,n) (—m,—n) and complex conjugation. Consequently, we need only numerically solve 
for the m > modes — although we still have to solve for both positive and negative n modes (except for m — 0). The 
computational burden is further reduced by noting that for I + m = odd the spheroidal harmonics vanish at 8 = n/2 
and as a result the corresponding ip^ , as well as their t, r, (p derivatives, automatically vanish at the particle (the 6 
derivative is also zero there, from symmetry). 



A. Numerical boundary conditions for the homogeneous fields 

The asymptotic boundary conditions for the radial equation are presented in Eqs. (j26[) and (|27[). However, in order 
to solve the radial equation numerically, boundary conditions need to be specified at finite radii. We denote the r* 
radius of the inner boundary by r* ln <C — M and that of the outer boundary by r* out 3> M (we discuss how these 
radii are chosen in practice in the algorithm section below). In order to determine the boundary conditions at r„ in 
and r* ut we use the ansatze 



k, 



= e+ "* out E4^;, (52) 

fc=0 

^ ^ = e- w ""£c fc -(r in -r + ) fc , (53) 



fc=0 



where ri n — r(r*; n ), J"out — r ( r *out), and the truncation parameters fci n , ut are chosen such that the boundary conditions 
reach a prescribed accuracy (see algorithm section below) . Substituting the above series into the radial equation (|24)) 



gives recursion relations for the c^ >0 coefficients in terms of Cq , respectively. In practice we take = 1. The explicit 
form of the recursion relations is rather unwieldy and can be found in appendix C of Paper I. 

In the case of the static m = n = modes there exist simple analytic solutions for ip^ , determined by regularity 
at the event horizon and at spatial infinity; they read 

^ho = a - rP i( z ) > tfoo = a + r Qi( z 1 ' ( 54 ) 

where a_, a+ are arbitrary constants, and are the Legendre function of the first and second kind respectively, 
and 



M 2 + a 2 

, a a (r-M). (55) 



B. Algorithm 

We now outline the explicit steps in our numerical calculation. 

• Orbital parameters. For given spin parameter a, orbital eccentricity e and semi-latus rectum p, calculate all the 
necessary kinematical entities associated with the geodesic orbit (£, £, fi r , T r , etc.) using the formulae given 
in section [TTJ 

• Boundary conditions. For a given Imn mode calculate the boundary conditions using Eqs. (|52[) and (|53[) with 
Cq = 1. For both boundaries we choose & max such that the relative magnitude of the fc max + 1 term drops 
below a given threshold which we take to be 10 -12 . In order to solve the radial equation (|24|) we first need to 
numerically invert Eq. (1231) to get r(r*). Machine accuracy places a limit on the smallest r* that the inverter 
can resolve and we take this value to be the inner boundary of our numerical domain r*j n . This value depends 
on a. For a = we find r*j n ~ — 65M and for a = ±0.998Af we find r»; n ~ — 262M. 

There is a fair amount of freedom in the choice of the location of the outer boundary r» out . Calculating the 
boundary condition using Eq. (|52l) is computationally cheap in comparison to solving the radial equation, and 
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therefore it is worth trying to place the outer boundary as close in as possible. In practice, for each Imn we 
initially attempt to compute the boundary condition with r wnt = 1000A/. For most modes this is adequate 
but for a few modes we find that the terms in the series start to grow after initially decreasing. If we detect 
this behavior we move r* ou t out by 5000ikf and try again. We repeat this procedure (at each stage moving the 
boundary out by 5000M) until the series converges to within the required accuracy. 

Homogeneous solutions. Using the boundary conditions as determined above we numerically solve the radial 
equation (|24[) for ipj 1 using the Runge-Kutta Prince-Dormand (8,9) [gsl_odeiv_stepjrk8pd] method from the 

Gnu Scientific Library (GSL) [29|. As we are using the method of EHS we solve for the outer field ibf between 

1 1 lmuj 

r* ut and r* m i n and for the inner field ib7 from r*- m to r* max . We store the value of the fields ib^ and their 
Ttr derivatives at a dense sample of points between r m ; n and r max , equally spaced in x- 

Inhomogeneous solutions. The next step is to compute the scaling coefficients using Eq. (|45|) . In prac- 

tice we calculate the integral in ([4"5"j) in terms of \ rather than t, using Eq. (|A2[) of Appendix [A] to con- 
vert between the two. The integration is performed numerically using the standard QAG adaptive integrator 
[gsl_integration_qag] (with the 61 point Gauss-Kronrod rule) from the GSL. The integrator automatically 
requests the values of the integrand between r m ; n and r max (or equivalently between x = and X — 7r ) required 
to perform the integral to within a set relative accuracy of 10~ 12 . The values of ip^ requested by the integrator 
are generated by locating the nearest of our sampling points stored between r m ; n and r max and using the value 
of the field and its derivative stored in the previous step as input to the Runge-Kutta algorithm in order to 
calculate the value of the field at the requested point. 

Determining n max - Although the sum in Eq. (|49p is formally over all n, in practice we of course sum over a 
finite number of modes, with |n| < n max . In our code we calculate the n = mode, then the n = —1, 1, —2, 2 . . . 
modes (in this order), until the relative n-mode contributions to both the field and its r derivative drop below 
a given threshold (which we take to be 10 -11 ). Typically we find that the contribution from the negative n 
modes drops below the threshold before the equivalent positive n mode does, especially for higher eccentricities 
(a similar behavior has been observed by Hopper and Evans [15j]). 

Spherical-harmonic projection. Once we have computed all the required Imn modes up to some I = /max (see 
below) we construct the extended spherical-harmonic <^>t modes using Eq. (|49p . The actual time-domain Z-mode 
is then constructed using Eq. (|51l) . and the l-mode contribution to the full force is obtained via Eq. (j52")l . When 
(i^O the coupling between the spheroidal and spherical harmonics implies that some of the spherical-harmonic 
modes with / < i max will have significant contributions from uncomputed spheroidal modes with I > Z max . We 
denote by Z max the highest spherical-harmonic mode for which the uncomputed spheroidal modes have a total 
relative contribution < 10~ 12 . 

Large-l tail contribution. To calculate the t, r, <p components of the SSF it is convenient to split them into 
conservative and dissipative pieces using Eqs. (|4"P|) . Regularization of the two pieces is then done using Eqs. ([57)1 
and p8p. The dissipative component requires no regularization and it converges exponentially fast, whereas 
the conservative piece converges like cx Z _1 — see Fig. Q]for an example. For a typical Z max = 25 the dissipative 
component is computed to a high degree of accuracy, but the slow convergence of the conservative piece neces- 
sitates an extrapolation of the regularized modes to account for the uncalculated large-Z tail. Our method for 
calculating the large-/ contribution is presented in section IV. C. of Paper I and we implement it here in the 
same form. We typically find the contribution from the first few I modes to be opposite in sign compared to 
that of the remaining modes, which sometimes results in that the contribution from the extrapolated / > 25 tail 
becomes as large as that of the calculated I < 25 part of the sum. The error from extrapolating the high-Z tail 
of the mode sum is by far the largest source of numerical error in the calculation of the conservative SSF. 

Once the conservative and dissipative pieces of the SSF have been obtained, the full SSF is calculated by simply 
adding the two together. Constructing the full SSF this way (i.e., conservative and dissipative pieces in separate) 
has practical advantages. For example, the behavior of the large-Z contributions to the full component F t near 
the orbital turning points transitions from an exponential decay (at smaller V) to an oc l~ 2 fall off (at larger I), 
as demonstrated in Fig. [2J It is then difficult to estimate the tail contribution to the full SSF component F t 
directly. However, when separated, the conservative and dissipative pieces of F t exhibit clear power-law and 
exponential behaviors, respectively. 
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Figure 1: Behavior of the Z-mode contributions to the mode sum for Ft (left panel) and F T (right panel), shown separately for 
the conservative and dissipative pieces, as well as for their sum. Here we have set (a,p,e) = (0.5M, 10M, 0.2) and the SSF is 
calculated at \ — t/2. The straight solid line is a reference line oc l~ 2 . From standard mode-sum theory we expect the Z-mode 
contributions to the full SSF and to the conservative component of the SSF to fall off as l~ 2 for large I. The curved solid 
line is an exponential reference line. We expect the i-mode contributions to the dissipative component of the SSF to decay 
exponentially with I. Numerical round-off error dominates the behavior of the dissipative modes for I > 15 
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Figure 2: Same as in Figure [T] this time showing the i-mode contributions to the SSF at \ = 0.010472 (= 27r/300), i.e., 
very close to periastron at r — r m i n . Near the orbital turning points the contribution to the full Ft (triangles, left panel) 
transitions from an exponential decay to an oc l~ 2 fall off. This transition (in this example around I ~ 8) makes extrapolating 
the large-Z tail of the full SSF difficult in this portion of the orbit. This problem is completely avoided if we separate the SSF 
into conservative and dissipative pieces, calculate the large-Z tail contribution to the conservative piece, and then add the two 
pieces together to recover the full SSF. Similar behaviour is observed for the F v component (not shown). No transition is 
observed for F r (right panel) as in this case the conservative piece dominates the SSF even near the orbital turning points. 



VI. CODE VALIDATION AND RESULTS 



The analysis of large-Z behavior, demonstrated above in Figs. [T] and [5J gives an important internal validation test of 
our numerical procedure and of the code itself. This is particularly so for the conservative piece, where the observed 
l~ 2 fall-off of the modal contributions rely on a delicate cancellation of as many as 3 terms in the l/l expansion [the 
0(1), 0(1) and 0(Z _1 ) terms]. This test obviously probes only the large-Z contribution to the SSF. A complementary 
test, which probes primarily the low-Z portion of the mode sum, is achieved by comparing the work done by the 
(dissipative piece of the) numerically computed SSF with the asymptotic fluxes of energy and angular momentum in 
the scalar radiation, which we may extract independently from the numerical data. 

In this section we present sample results from our SSF code, and in particular discuss the above flux comparison, 
which we consider a strong quantitative validation test. Since our work represents a first full implementation of the 
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method of EHS in a SF calculation, we then take some time (in subsec. IVI C|) to examine the efficiency of our method. 

A. Dissipative piece of the SSF 

In Paper I we derived formulae for the average flux of scalar energy radiated to null infinity and through the event 
horizon. The formulae we derived specialized to circular equatorial orbits but the method of derivation was quite 
general and is readily extended to eccentric equatorial orbits. The resulting equations are given by 

(E + ) = ^E^I^J 2 ' = ^E^«(^n-m0 + )|(7r n p , (56) 

Imn Imn 

( L +) = T-Y.muimnlC^ | 2 , (£-) = ^— VmK m -mfi + )|C: | 2 , (57) 

v ' 47T ^ — ' Imn' ZTTTa- ^ — ' Imn 1 v ! 

Imn Imn 

where the + and — subscripts denote the flux radiated to infinity and down the event horizon, respectively, f2+ = 
a/(2Mr + ), an overdot denotes differentiation with respect to t, and (•) denotes a t-average over an orbital period T r . 
The amplitude coefficients are given by 

Cf =Cf 4, (58) 

Imn Imn u 

where, recall, Cp are given in Eq. (pPfj) . and are the leading coefficients from Eqs. (I52"j) and (I5"3"l) . From the form 

of (E_) and it can be seen that the Zmn-mode contribution to the horizon flux becomes negative whenever 

Umn < mQ + . For these superradiant modes (30j the scalar particle gains energy and angular momentum at the 
expense of the mass and rotation of the central black hole. 

Using Eqs. ([7]) and ([5]), and noting that the conservative contribution to the averages {£) and (£) vanishes by virtue 
of Eq. (|40|) . we obtain the relations 

-^WAmJ - (59) 

where we have used Ut = —£ and u v = C, and where A/i is the net change in the particle's rest mass over a period 
T r . The latter is found from Eq. © to be identically zero, 

A M = 0, (61) 

since in our case 4> fi (r) comes back to itself after a period T r . The orbital energy and angular momentum dissipated 
by the SSF over a period T r should be balanced by the total energy and angular momentum radiated to infinity and 
down through the event horizon over that same period, i.e., 

-H(S) = (£) total = (£+) + (£_), (62) 
-H{£) = (L) total = (L+) + (L-). (63) 

We used our code to calculate both sides of the balance equations (1521 and (|63l) for a variety of orbits and black holes 
spins. Table U summarizes a small sample of our results. In all cases considered we found a good agreement between 
the local dissipative SSF and the radiated fluxes. This comparison tests primarily the low-/ modal contributions to 
the SSF, because the amplitude of these contributions falls exponentially with I. 

B. Sample results 




Using the algorithm outlined in section fVBI we have calculated the SSF for a variety of black hole spins and orbital 
parameters. The highest eccentricity we have been able to explore is around e = 0.7 (see below for a discussion of 
limiting factors). In Fig. [3] we present sample results for the SSF along a variety of orbits and in Fig. Uwe show an 



14 



a/M p/M 


e 


(£) x »(M/q) 2 1 


-H{£)/(£)total| 


{£) X fiM/q 2 1 


- M |(>C)/(i}total| 


0.9 


10 


0.2 


2.6862422 x 10~ 5 


-7.4 x 10" 8 


-8.3593539 x 10 -4 


-7.4 x 10" 8 


0.9 


10 


0.5 


2.4856622 x 10~ 5 


9.5 x 10" 8 


-6.2962019 x 10" 4 


7.2 x 10" 8 





10 


0.2 


3.213314062 x 10~ 5 


1.2 x 10~ 10 


-9.62608845 x 10~ 4 


1.0 x 10~ 10 





10 


0.5 


3.329332 x 10~ 5 


1.6 x 10" 7 


-7.844684 x 10~ 4 


3.3 x 10~ 7 


-0.5 


10 


0.2 


3.65656098 x 10~ 5 


2.8 x 10~ 10 


-1.06932319 x 10~ 3 


2.0 x 10~ 10 


-0.5 


10 


0.5 


4.33567 x 10"" 5 


3.1 x 10" 6 


-9.70033 x 10~ 4 


2.4 x 10" 6 


0.2 


6.15 


0.4 


3.42797 x 10~ 4 


2.5 x 10 -6 


-3.92668 x 10" 3 


2.2 x 10 -6 



Table I: Orbital energy and angular momentum dissipated by the SSF and comparison with the radiated fluxes, for a variety 
of orbits with p = 10M. The last row shows data for a "zoom- whirl" -type orbit (cf. Fig. [4}. The average dissipation rates {£) 
and (£) (4th and 6th columns) are calculated from the local SSF using Eqs. (|59[1 an d (|60[) . The radiated energy and angular 
momentum (-E) to tai and (L)totai are extracted independently from the asymptotic fluxes using Eqs. (156[) and (|57[) . The relative 
differences displayed in the 5th and the last columns verify that the balance relations (I62|l and (I63|l are satisfied. We believe 
the dominant source of residual discrepancy comes from the numerical integration in Eqs. (I59|l and (160[) . 



example of a 'zoom-whirl'-type orbit. Table [TTI displays some numerical results for the SSF for various orbits and spin 
values. 

In the Schwarzschild case (a — 0) it is possible to compare our results with those from the recent analysis of 
Canizeras et al. |36| , who used a pseudospectral algorithm formulated in the time domain. We find a good agreement — 
see Appendix [B] We have also tested the output of our code (in the a = case) against more detailed (unpublished) 
data from a time-domain code by Haas (35j | . 
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Table II: Numerical results for the dissipative and conservative pieces of the SSF for a sample of orbits. The full SSF is obtained 
by adding the two pieces together. The SSF is sampled at \ values corresponding to the points marked along the orbits in Figs. 
[3] and [4] Data for the conservative components include an estimate of the uncertainty from the large-/ extrapolation, which 
dominates the overall numerical error in these components; this is indicated by figures in brackets, showing the uncertainty 
in the last quoted decimal. We used the method described in Paper I to estimate this error. In the case of F^ iss no large-Z 
extrapolation is needed, and the accuracy is much improved; in this case we believe all figures shown are significant. The SSF 
data for this table was obtained with typical values of i max between 15 and 20. 



C. Computational performance 

The computational burden increases moderately with \a\ and more rapidly with e. The larger the spin magnitude 
\a\ is, the stronger the coupling between spheroidal and spherical modes becomes, and the more I modes need be 
calculated for given ' max . The higher the eccentricity e, the broader the Fourier spectrum becomes and the more n 
modes need be calculated for each l,m. Moreover, larger eccentricity also leads to a stronger spheroidal-spherical 
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Figure 3: Sample SSF results. The top left panel shows orbits with (p, e) = (lOAf, 0.5) in the equatorial plane, for three 
different black hole spins: a = (dotted, blue curve), a = 0.9M (solid, red curve) and a = — 0.5M (dashed, green curved). For 
each orbit we show one complete revolution, from one periastron to the next, with markers indicating the points taken for the 
sample data in Table [Til The other three panels show (reading clockwise) the F r , F v and Ft components of the SSF, for the 
three orbits shown in the top left panel. \ = corresponds to a periastron passage. 



r sin(<^) /M 




Figure 4: Shown in the left panel is a 'zoom- whirl'- type orbit with parameters (a,p,e) = (0.2M, 6.15M, 0.4). The right panel 
shows the corresponding components of the SSF along the orbit, x = is a periastron of the orbit. Markers indicate the 
location of the data points shown in Table [TT| 
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Figure 5: Computational cost. We show the total time required to calculate up to (a fiducial) / max = 15, leading to fractional 
accuracies of order ~ 10~ 4 in the SSF. In the Schwarzschild case (a = 0) we can compute the SSF for eccentricities up to e = 0.7 
in around 12 hours. For a = 0.9M the calculation requires more time, primarily because the coupling between the spheroidal and 
spherical harmonic modes necessitates calculation of higher spheroidal harmonic modes, which is computationally expensive. 
At low eccentricities, our frequency-domain algorithm is vastly faster than any existing time-domain method. 



coupling, because the spheroidicity parameter a 2 that determines the strength of this coupling is proportional to 
w mn' which is larger for higher n harmonics. Using the current version of our code we were able to explore spin 
parameters in the range — 0.99Af < a < 0.99M and eccentricities in the range < e < 0.7. Beyond these ranges the 
computational burden becomes prohibitive. 

In Fig. [5] we plot the CPU time required to compute the SSF on a standard desktop machine (dual-core, 3GHz). 
We used a fiducial i max = 15, giving SSF fractional accuracies of order ~ 10~ 4 . We show results for a = 0.9M and, 
for comparison, a = 0; results for a = — 0.9M are found to be similar to those for a = 0.9M. Note that the data 
for a = probes the performance of the EHS method, while the Kerr results also reflect the increased computational 
burden due to mode coupling. 

In the Schwarzschild case (a = 0) we find that for e < 0.4 the computation time grows only ~ linearly with e, 
somewhat more rapidly at higher eccentricities, and very fast for e > 0.6. Nonetheless, the required computation time 
for e = 0.7 is still only around 12 hours. For comparison, an equivalent time-domain computation [3l| (on a similar 
machine and with similar accuracy standards) takes several days. At lower eccentricities our gain in speed/accuracy 
is very substantial. Our results for a = highlight the efficacy of the EHS method. 

In the Kerr case mode coupling adds to the computational burden, and more so with growing eccentricity. As an 
example, for a = 0.9 with e = 0.2, we find that the spheroidal-harmonic I — 15 mode has significant contributions 
from all tensorial-harmonic modes 8 < I < 22. This results in a more rapid growth in CPU time as a function of e, 
compared to the Schwarzschild case. For spins as high as \a\ = 0.9, eccentricities greater than ~ 0.5 are practically 
beyond reach for our current code. However, for small eccentricities our algorithm is extremely efficient even at high 
spin. 



VII. ISCO SHIFT 



The ISCO shift due to the conservative piece of the SSF for a particle in orbit about a Schwarzschild black hole 
was first calculated by Diaz-Rivera et al. [19|]. More recently the ISCO shift due to the conservative piece of the 
gravitational SF for a similar orbital setup has also been calculated [l|, [l6[ . Here for the first time we calculate the 
conservative SSF correction to the innermost stable equatorial orbit (ISCEO) for a particle in orbit about a Kerr 
black hole. The following derivation follows closely that of Ref. [l|, but we adapt it here to Kerr spacetime. 
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A. Test particle case 

In this section we review the notion of the ISCEO for a test particle in Kerr geometry, ignoring all SSF effects. We 
begin with the radial geodesic equation for equatorial orbits in Kerr spacetime, given by [26j 

= ^ { W + a2 ) - «4 2 " A K + - = nr P M . (64) 

Differentiating with respect to r gives 

d 2 r 1 dIZ 

— s -=J r e s{r p ,£,£), Fes{r p ,£,£) = -— , (65) 
dr z 2 ar p 

with J^ff being an effective radial acceleration. 

For slightly eccentric geodesies (e C 1) we can expand the particle's radius as a function of r in the form 

r-p(r) = r + en(T) + 0(e 2 ) , (66) 

where ri(r) is independent of e and comparison with Eq. (|S|) allows us to identify ro = p. In this equation and 
throughout this section we use a subscript '0' to denote the circular-orbit value (e = 0) and a subscript '1' to denote 
the 0(e) perturbation in the quantity's value, holding ro fixed. Substituting Eq. (|66[) into Eq. (|65p and reading off 
the 0(e) terms we find 



d 2 r 1 _ dF cS {r pi £,C) 



dr 2 dr p 



n , (67) 



e=Q 



where we note £ \ = d = by virtue of the quadratic dependence of £ and C on e [recall Eqs. ([3]) and dl]), replacing 
p with tq]. Hence the 0(e) radial motion is a simple harmonic oscillator, 



with frequency 



where we have denoted 



dT cS (r p ,£,C) 



dr p 



= -u r n , (68) 



M[r {r -6M) + 8avr -3a 2 } 



e=0 



r4(r - 3M + 2av) 



v = y/M/ro ■ (70) 

Assuming a periapsis passage at r = we obtain n = — ro cos(oj r r) and hence 

r p (r) = r (l - e cos w r r) + 0{e 2 ) . (71) 

The location of the ISCEO is defined to be the radius r; s for which w r (ri s ) = 0. Solving the quartic in Eq. (|69l) . 
and defining a = a/M, one finds 

r is = Af{3 + K-sign(a)[(3- 7 )(3 + 7 + 2«)] 1 / 2 } , (72) 

where 

7 = 1 + (1 - a 2 ) 1/3 [(l + a) 1/3 + (1 - 5) 1/3 )] , (73) 
k = (35 2 + 7 2 ) 1/2 , (74) 

which to the best of our knowledge first appeared in Ref. [25j]. For the Schwarzschild and the extremal (|a| = M) 
prograde/retrograde cases the unperturbed ISCEO is located at 6M, 1M and 9M respectively. The azimuthal orbital 
frequency of a test particle at the ISCEO is given by [26] 

M 1 / 2 

^ - r 3/ 2 + aM y 2 ■ (75) 

Recall our convention Q v > 0, with a > for prograde orbits and a < for retrograde orbits (this is different from 
the convention of Ref. [26| ) . 
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B. SSF correction to the ISCEO 



Let us now derive the conservative SSF correction to the ISCEO location. The perturbed equations of motion 
including conservative-only SSF effects are 



df 



1 p i (cons) 
'H 1 r t 



dC 

df 



df 2 



F eS (f p ,e, C)+fi- 1 F 



l^i_L(cons) 



_L(cons) ' 



(76) 
(77) 



where hereafter we denote perturbed quantities by an overbar, and we define £ = —u t and C = u v (no longer 
necessarily conserved along the orbit). We use the sub/superscript _L(cons) to denote the conservative piece of the 
SSF perpendicular to the particle's 4- velocity [see Eq. ((7} and section [Til A] . 

We assume that the radius f p (t) of the SSF-perturbed slightly-eccentric orbit can again be formally expanded about 
a circular orbit of radius ro, 



f p (r) = r + efi(T) + O(e 2 ) , 
where fi depends of ro but not on e. We similarly expand 

£ = £ + eSi(r) + 0(e 2 ), C = C + eA(r) + 0(e 2 ), 



(78) 



(79) 



where £q and Co are the SSF-perturbed values of £q and Co along the circular orbit of radius ro. To find Co and £q 
we simultaneously solve df/df = and d 2 f/d? 2 — [hence lZ(f p ,£ ,£, ) = with d1Z/df p (f p , Bq,Cq) — 0]. This gives 



£ a = (l-'Sv 2 + 2dv z )- 1 ' 2 
Co = r {l-3v 2 + 2dv 3 )~ 1/2 



1 - 2v 2 + dv 3 - —F 



2fi 



±o 



v(l - 2dv A + d z v^) - 



r (l 



2vfi 



-FT 



(80) 
(81) 



where Fj_ is the circular-orbit value of fj\ co > (note the r component of the SSF is purely conservative along a 



circular orbit, so the label 'cons' becomes redundant in this case). 
The 0(e) part of Eq. (f77f now takes the form 



d 2 n 

dr 2 



-u} 2 n , 



where 



LO r . = — - 



dr„ 



T eS (f p ,£,C) + ^ 1 F'l 



_L(cons) 



(82) 



(83) 



Here £ , C and F T , 



±(cons) are thought 01 as functions of f p along the orbit (for given ro,e), and the f p derivative is 
taken with fixed ro, e. The form (|83|) assumes that £ , C and -FX (cons) depend explicitly on e (when ro and f p are held 
fixed) only through e 2 . That this is true in the Schwarzschild case was shown in Ref. (l6j based on a simple symmetry 
argument, and the same argument carries over to the Kerr case. The perturbed radius is thus a simple harmonic 
oscillator in r with frequency ui r (for Co 2 > 0) , and choosing t = at periastron passage we have 



r, p (t) = ro (1 — ecosw r r) . 
Using fi = — ro cosw r f and recalling Eq. (|76j) . we may now write 

FI (co „ s) = Flo + e*Ii cos Q r f + 0(e 2 ) 



ci_L(cons) 
77,_L(cons) 



euJrF^ 1 sino; r T + 0(e 2 ) , 
euJrF^- 1 sinw r f + 0{e 2 ) , 



where we have defined 



Fli = -ro 



"•^I(coas) 

dr„ 



dC 

df„ 



F 



±i 



d£ 

dfn 



(84) 

(85) 
(86) 
(87) 
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Then, using these definitions in Eq. (|83l) and substituting for 8q and Cq from Eqs. (|80| and (f8TJ) . we obtain 
-2 2,*L 3 - 12i- 2 + 9^ 3 F[ aM + a 2 v + ro(r -3M)v F^ 

U r = W r + — — 7i o„.2 I o^..^— " 2 Q/2 i == = ~ t 89 ) 



-2a 



r /J r (l - 3w 2 + 2av 3 ) (j, r^ /2 ^/r - 3M + 2av M 

a(r + M) + a 2 u - 3Mr v F^ 1 



— , .2 , A, .2 



i"o /2 vVo - 3M + 2av V 



uj r + Au} r (r ,a) , 

where Aw 2 (ro,a) denotes the terms proportional to the SSF and we have neglected terms quadratic in the SSF. 
Near the ISCEO, the unperturbed frequency (squared) may be expanded in the form 

u 2 {r , a) = A{a){r - r is ) + O(r - r is ) 2 , (90) 

where, recall, r; s denotes the location of the unperturbed ISCEO, given in Eq. (1721) . and 



., . du: 2 
A{a) = — 1 
or 



= 3M rj + M 2 r is (18 - 55 2 - 385^) - 7a 2 M 3 (~av is ~ 4) + Mrf s (irav is - 10) 
r =r la '* rf s + M{2av is - 3)] 2 



with Vi s = v M/vis- By definition, w 2 vanishes at the (shifted) location of the ISCEO: w 2 (ro = ris) = 0. By 
substituting Eq. (|90|) into Eq. (|89|) . setting r = f; s and u) r = 0, and solving for fj s at linear order in the SSF, we find 
the SSF-induced shift in ISCEO radius to be [through 0(q 2 )} 

An s = f is -r is = -^4^. (92) 
A{a) 

Note on the right-hand side of Eq. we have substituted r; s for f- ls as this term is already of 0(q 2 ). When a = 0, Eq. 
(1921) reduces (upon replacing the SSF components with the gravitational SF components) to the ISCO shift formula 
found in Ref. 1], namely 

Ar is (a = 0) = (M 2 /^) (216Fl 0is - 108*I lta + V3M~ 2 F^ S ) , (93) 

where the 'is' subscript denotes a quantity's value at the unperturbed ISCO. 

We next consider the conservative SSF shift in the azimuthal frequency The frequency associated with the perturbed 
circular orbit of radius r — rp is given by 



difip _ dip p /d,T _ g$ v £ - gjf £ 
dt ~ di/df ~ g^Co-g^So 



where are the background metric functions evaluated on the perturbed circular orbit. Substituting for Eq and Cq 
from Eqs. (|80|) and (|8TT) . taking r*o = r; s + Ari s and keeping only terms through 0(q 2 ), we find the relative frequency 
shift at the ISCO to be given by 



Afi vis _ fi v i s - 3Ar is rf s (n s -3M + 2avi s )fi 1 F r 



_L0is 



^t^is ^ipis 2(7*is + av is ) 2M(n s + av is )(rf s - 2Mr- lt 



(95) 



For a = the above formula reduces (when the SSF components are replaced by the gravitational SF components) 
to that found in Refs. 0, [lj|, namely 



AR„ is , Ar ifi 27M 



^( a = 0) = -=^-^l J F[ n . . (96) 

O vis 1 7 4M 2/j J - 0ls v 7 

All that remains is to rewrite the above expressions for Ar- ls and A£l V i S in terms of the full Boyer-Lindquist 
components of the SSF (rather than the normal components F^)- Specifically, recalling Eqs. and (|95p . we will 
need expressions for F^ ± , Fj_ x , F^ 1 and F^ 1 in terms of the quantities F®, F^, F^ and F$ arising, in analogy with 
Eqs. (|85 |l -(|87 p . from the formal e-expansion of the full conservative SSF: 

F™ ns = F° + eFtcosuj r T + 0{e 2 ) , (97) 
F™ ns = eu, r Fismu, r T + 0(e 2 ) , (98) 
F™ ns = eu r Flwa.u r f + 0(e 2 ) . (99) 
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(We prefer to work with the covariant components of the SSF, which are the ones returned by our code.) Starting 
with the radial component, we write 



Fl (cons) = (9 rp +u r u^F^ s , (100) 

and formally expand both sides of the equation in e using Eqs. ([H5"]l and (l^7 )) - ([M)) . noticing u r u^Fp° ns = 0(e 2 ). 
Comparing the 0(e ) and 0(e l ) terms on either side then readily gives 



Flo = 9l r F° r , (101) 
F r xl = So" r 



F r ° . (102) 

r=r 



dr 

For the t and ip components we similarly obtain 

F^ 1 - (1 + Clg^ - £ C g^)Fl + (C 2 ^ - £oCo9$)F? + C r F? , (103) 
F t xi = (1 + - E Q C^)F} + {E 2 g^ - S C g^)F^ - £ r Q F? , (104) 

where we have also used u r — erQUj r smuj r f [recall Eq. (|84|)] . 

The shifts in the location and frequency of the ISCEO can now be constructed from the full SSF by substituting 
Eqs. (|101l) - (|104j) (evaluated at the ISCEO) into Eqs. (j§2"]) and (j5B"j) . The resulting formulae are cumbersome so we 
leave them implicit. For a — the formula for the radial ISCEO shift is found to reduce to 

Ar is (a = 0) = (M 2 / M ) (216F° ls - 72F^ + 6V2F t \ s + ^^ is ) . (105) 

which is in agreement with Eqs. (51) of Diaz-Rivera et al. (l9| . 



C. Numerical results 



In order to implement Eqs. f[92]) and (|95"T) we require the values of F®, F*, F* and Fh, all evaluated at r = r is . 
The first piece of data, F®, is simply the radial SSF component evaluated along a circular equatorial orbit of radius 
fo = Ti s , and we can compute it with great accuracy using the circular-orbit code of Paper I. The computation of the 
other quantities, which are associated with a slightly eccentric orbit, is more delicate. Recalling Eqs. (|97|) -(|99 l) we see 
that they may be extracted using 

F^(p,e) = 2uj r (eTT)- 1 F r cons cos(o; r r)dr , (106) 

Jo 

Ffo,e) = 2(e7r)- 1 f F- I1S sin ( Wr r)dr , (107) 
Jo 

where a £ {t,w} [we are allowed here to remove the bars off ui r and f since the quantities F a are already 0(/i)]. As 
noted in Ref. |l|, both limits can be taken simultaneously by picking points along a suitable curve in the e-p plane. 
As also discussed in Ref. [1], for our 0(e)-expansions to be valid we must have both e<Cl and e <C (p— r- ls )/M along 
the curve (so that we keep sufficient distance from the separatrix as we approach the ISCEO). The final result should 
be independent of the particular path taken through the e-p plane and we use this fact as a validation test of our 
numerical scheme. In practice we calculate F[ , F^ and F* at various points along three curves given by p = ri a +M\fe, 
p = fi s + %M\fe and p — r ls + Me 1 / 3 , and then extrapolate each set of data to the ISCEO — see Figs. [6] and [3 The 
(small) variance in the extrapolated values from the 3 curves is used as an error estimator for the iq's. 

Once the Fi's are at hand, we use Eqs. (|§2"j) and (jHU) to compute Arj s and AJ7 v i s for a variety of a values. The main 
source of error in our final results comes from the e — > extrapolation involved in extracting the F\ functions (the 
error in Fq 1s is relatively much smaller and can be neglected). As a rough estimator of the error in An s and Att^s 
we use the variance in the values of these quantities when using the three different extrapolation curves mentioned 
above. [Notice we do not consider here the more conservative estimate obtained by adding up the contributions to the 
error from the various -Fi's (either in absolute values or in quadrature) since these are clearly correlated; we believe 
our estimate better represents the actual uncertainty in the final results.] 

Our results are presented in Table IIIII and Fig. [5] We observe that Arj s increases monotonically as the black hole 
spin is varied from a = —0.9M to a = 0.9M , and it changes sign from negative to positive around a = 0.8Af . The 



FX = hm ]im Fr(p,e) , 
^iis = lim limF^(p,e) , 
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Figure 6: Calculation of F^ is , Fjr ls and Fh s by extrapolation along three paths in the e-p plane for the case of a Schwarzschild 
black hole (a = 0). The sampling paths in the e-p plane are shown in the bottom right panel. The other three panels show 
numerical data points for F^, and Ft as extracted from the conservative piece of the SSF using Eqs. (|106l) and (| 10T[I . Solid 
curves are cubic interpolations of the numerical data points, and the extrapolated values at p = 6 (which, in theory, should 
not depend on the choice of curve) represent our numerical predictions for F^ is , F^ and F^ ia . The small variance in these 
extrapolated values serves as a rough measure of error. The thick dot on the vertical axis marks the values found by Diaz- Rivera 
et al. 19|. For F* ia and F^ iB we find a close agreement with their results, but for _F t lls there is a discrepancy at a level (~ 2%) 
which we cannot explain. (The code used by Diaz-Rivera et al. cannot be retrieved to allow a careful examination of this 
discrepancy [32]]: we are, however, quite confident in our results given the the good agreement between the values extrapolated 
from the different curves.) 



relative shift in the azimuthal frequency at the ISCEO, AO^s/fi^is, is always positive between a = — 0.9M and 
a = 0.9M. For retrograde orbits the relative frequency shift remains similar to that found in the Schwarzschild case 
(a = 0), while for prograde orbits it decreases rapidly with increasing spin a. Our code is not sufficiently accurate to 
explore the near-extremal case, so the behavior of Af^is (and of Ar- ls ) there remains unclear. 



VIII. VARIATION OF REST MASS 



As discussed in Sec. Ill Al the SSF has a component tangential to the particle's worldline, which leads to the particle 
having a dynamically varying rest mass (20| . (This situation is special to our particular SSF theory; in the equivalent 
electromagnetic and gravitational cases the rest mass is conserved. It is possible to construct a scalar field theory 
where the particle's rest mass is conserved but only at the cost of making the field equation non- linear (20|.) Previous 
studies of this phenomenon in cosmological spacetimes [33l . |34| have found a range of possibilities including a periodic 
mass variation as well as cases where the mass dissipates entirely. 

In our setup, where the motion is intrinsically periodic, the field returns to its original value after one orbital 
revolution and thus from Eq. we see that the net change in the particle's rest mass will be zero. Furthermore, 
examining Eq. ([8} and recalling the symmetry relations expressed in Eq. (|40p . we can see that d/i/dr is symmetric 



22 



-0.015 I r- 

-0.02 

-0.025 - 
-0.03- 

^-0.035 - 

X -0.04 - £ 

<ClH -0.045 - 
-0.05 

-0.055 - 
-0.06 



2.3 




2.4 



2.5 



V = r-m + -JeM o 
p = r i3 + e 1/3 M □ 
p = r is + 
~2£ 27 Z8 23 
p/M 





0.125 - 




0.12 




0.115 




0.11 - 












0.105 


«l 




X 


0.1 - 


<^ 






u.uyo 




0.09 - 




0.085 - 




0.08 - 



2.3 



V : 



r ls + VeM 
p = r is + e 1/3 M 
p = r is + 3/2VeM 




2.4 



2.5 2.6 

p/M 



2.7 



2.8 



Figure 7: Same as in Fig. [6] but for a = 0.9M. We use the paths shown on the graphs to extract the values of F rls , F tls and 
F„ ls via extrapolation to the ISCEO, whose location at n s = 2.32088M is marked by the vertical line. (The results for Fh e , 
for brevity not shown here, are qualitatively similar to those of F^ is .) Sample numerical values for F^ is , Ftis an d -Pools > an( A the 
resulting ISCEO shifts for different spins a, can be found in Table uTIl 
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-0.1 


6.32289 


0.063294 


1.62329 x 10" 4 


7.98(2) x 10~ 4 


0.01833(2) 


-0.1141(2) 


-0.3 


6.94927 


0.055496 


1.31790 x 10" 4 


7.13(2) x 10" 4 


0.01504(5) 


-0.1050(5) 


-0.5 


7.55458 


0.049348 


1.27517 x 10" 4 


6.23(2) x 10~ 4 


0.01257(2) 


-0.0970(2) 


-0.7 


8.14297 


0.044372 


1.10762 x 10" 4 


5.36(2) x 10~ 4 


0.01068(3) 


-0.0909(6) 


-0.9 


8.71735 


0.040260 


9.60700 x 10 -5 


4.68(2) x 10~ 4 


0.00919(2) 


-0.08492(7) 



Ar i8 



fjM 



0.0309(5) 

0.004(1) 

-0.0226(7) 

-0.0471(2) 

-0.0651(6) 

-0.0804(6) 

-0.10070(4) 

-0.1172(2) 

-0.1268(1) 

-0.1350(3) 

-0.1505(9) 

-0.1620(4) 

-0.1723(8) 

-0.1802(4) 



-0.0107(3) 

0.0013(5) 

0.0118(7) 

0.01802(6) 

0.02203(5) 

0.0248(6) 

0.02779(1) 

0.02948(4) 

0.03020(3) 

0.03056(6) 

0.0309(4) 

0.03090(8) 

0.0306(2) 

0.02992(7) 



Table III: The conservative SSF effect upon the ISCEO location and frequency. Each row of the table corresponds to a particular 
value of the Kerr spin parameter a: the second and third columns show the values of the unperturbed ISCO radius n s and 
frequency fi^ia, and the fourth through seventh columns show the numerically-computed values of the SSF coefficients F^- iB , 
Frist ^tis an< i F<pis defined through the small-e expansion in Eqs. (|97[) - (|99[) . The last two columns display the SSF-induced 
shift in the radius and frequency of the ISCEO, as computed using Eqs. (|92[) and (|95[) . Figures in brackets are estimates of the 
numerical error in the last displayed decimals (in the data for _F° ia all figures are significant). Note the fractional error in the 
a — 0.8M results for Ari s and Af2 v i s is particularly large: this is a consequence of a delicate cancellation between the various 
terms in Eqs. Q92p and H96|) . which also leads to the vanishing of Arj s and Afi^ia at two (slightly different) spin values close to 
a = 0.8. For a = Diaz-Rivera et al. []J] obtained Ar is = -0.122701g 2 /^i and Att^/tt^ = 0.029l657q 2 /(fj,M). The small 
discrepancy is discussed briefly in the caption of Fig. [5] 



about the apastron and hence the rest-mass change from periastron to apastron (and visa versa) must also be zero. 
To within our numerical accuracy we observe this behavior in our data — see Figure O below. 

It is also interesting to examine how the rest mass varies along the orbit. The total rest mass change from periastron 
to a point with phase x along the orbit is given by 

MX) = -] F a SS (x)u a ( X ')^d X ' . (108) 

As illustrated in Fig. |H1 the particle's rest mass initially increases (though for zoom-whirl type orbits there is a slight 
decreases before the main increase) but then decreases so that the particle regains its original mass by the time it 
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Figure 8: The conservative SSF effect upon the ISCEO location and frequency. We plot here the numerical results shown in 
Table HITl as a function of the Kerr spin parameter a (curves are cubic interpolations). Vertical error bars indicate the estimated 
numerical error. Notice in the graph that two separate scales are being used for An s (left-hand scale) and Afi^is/f^is (right- 
hand scale). The radial shift is monotonically increasing with a and changes sign around a = 0.8M. The frequency shift 
similarly changes its sign (and turns negative) at large spin values. Note that although the change of sign in the radial and 
frequency shifts occur near the same spin value (a =~ 0.8M) the error bars on our results excludes the possiblity of the sign 
change occuring at the same spin value. The Schwarzschild ISCO shift results of Diaz-Rivera et al. 19] are marked (green, 
thick dots) for comparison. The small discrepancy is discussed briefly in the caption of Fig. [6] 

reaches apastron. Past the apastron the mass continues to decrease before increasing back to the original value at 
periastron. We also observe that the change in mass along the orbit is only weakly dependent on the black hole spin, 
for fixed (p, e). 

In our setup the particle's rest mass is conserved over an orbital period T r , but in a setup which allowed for the 
orbit to evolve through the action of the SSF this would no longer be the case (the field would no longer return to 
it's original value after one orbit). It would be an interesting project, which we do not pursue here, to consider the 
effect of the net mass loss on the insipral dynamics of the scalar charge. 

IX. SUMMARY AND FUTURE WORK 

In this work we have presented a first calculation of the SSF for a particle moving along an eccentric, equatorial 
geodesic of a Kerr black hole. Working in the FD we have made use of the recently proposed method of EHS to 
overcome the difficulties that would have otherwise hampered the convergence of the frequency mode-sum (in either 
Schwarzschild or Kerr spacetimes) . As this work represents the first complete SF calculation made using the method of 
EHS we have also explored the efficiency of the method and found it to perform extremely well. In the Schwarzschild 
case, calculating the SSF along an eccentric geodesic (at a fractional accuracy of < 10~ 4 , on a standard desktop 
computer) takes between a few minutes for small eccentricities and ~ 12 hours for e = 0.7 — still substantially faster 
than equivalent time-domain codes. Calculations in Kerr are more computationally demanding due to mode coupling 
(and in this case there are no time-domain codes at hand yet to allow comparison). 

In Table lU we provide sample SSF results, which may serve as comparison data for researchers working on alternative 
computation methods in Kerr. We performed several validation tests to establish confidence in our results. First, we 
checked that, for each point along the orbit, the l-mode contributions to the SSF exhibit the correct large-Z asymptotics, 
as predicted by mode-sum theory. Second, we compared the rate of energy and angular momentum dissipation by 
the SSF with the corresponding fluxes in the scalar waves radiated down the event horizon and to infinity. Finally, 
for the non- rotating case we compared our results with those of Haas [35[ and Cahizares and Sopuerta (36j. 
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Figure 9: Rest-mass change due to the SSF for a scalar charge in an eccentric equatorial orbit about a Kerr black hole. (Left 
panel) The variation in the rest mass, for orbits about a Schwarzschild black hole, as a function of x, is strongly dependent on 
the orbit's eccentricity. For e and p far away from the separatrix the mass initially increases and then returns to its original value 
at apastron before further decreasing and then once again returning to the original value at periastron. For a 'zoom-whirl'-type 
orbit the mass is observed to decrease slightly initially. (Right panel) Results for Kerr. The change in mass is weakly dependent 
on the black hole spin. 



We have used our numerical results to quantify two "physical" effects associated with the SSF in our scalar-field 
model. First, we calculated the shift in the ISCEO radius and frequency due to the conservative piece of the SSF, 
hence generalizing previous results by Diaz-Rivera et al. [HI to the Kerr case. Our main results for the ISCEO 
shift are summarized in Table Mil and plotted in Fig. El We then also examined, for the first time in a black hole 
spacetime, the variation in the rest mass of the scalar particle resulting from the component of the SSF tangent to 
the four- velocity. We confirmed that the particle exchanges mass/energy with the scalar field in such a way that (if 
the inspiral motion is ignored) no net mass loss occurs over a complete radial period. Sample results for the mass 
change are presented in Fig. [5] 

Several extensions of this work suggest themselves. First, one might generalize to generic, inclined orbits in Kerr 
geometry. We have made some progress calculating the SSF for circular inclined orbits (and will present our results 
elsewhere), but the challenge of generic orbits (ones that are inclined and eccentric) still lies ahead. Such generic 
orbits are tri-periodic, with the additional frequency coming from the longitudinal motion. This is certain to increase 
the computational burden and might further restrict the portion of parameter space that is tractable using our FD 
method. 

The important extension to the gravitational case is even more challenging. It is not clear if there is a second- 
rank tensorial generalization of the scalar spheroidal harmonics that would facilitate a full separation of variables in 
the frequency domain in Kerr geometry. In the absence of such formulation one could still tackle the gravitational 
perturbation equations in the frequency domain, by decomposing in the standard tensorial spherical harmonics and 
then properly accounting for the coupling between I modes that would result in Kerr. 

At least in the Schwarzschild case, however, it is clear that an EHS-based FD method is both viable and extremely 
computationally efficient. In this case, a standard tensor-harmonic decomposition can be employed to achieve a full 
separation of the gravitational perturbation equations (e.g., in Lorenz gauge), and one can proceed in a straightforward 
way to construct the gravitational SF using EHS in conjunction with mode-sum regularization, as in the scalar case. 
This method has already been implement recently by Akcay for circular orbits [ [ 9] , and work to generalize it to eccentric 
geodesies is well under way. 
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Appendix A: Some expressions entering the formulation of eccentric equatorial geodesies in Kerr spacetime 



We give here a few expressions that go into the calculation of t he g eodesies in section|TTJ The azimuthal angle f P (x) 
and Boyer-Lindquist time t p (x) along the geodesic are given by [22j 



Mx) = / d x — \~xi2, : 

Jo J(x',p,e)V r ' {x',P,e) 



t P (x) 



V t (x',P,e) 



dx'~ j .> 

o J(x',P,e)Vr' {x',P,e) 



where 



2Mx 2 

V r (x,P,e) = x 2 + a 2 + 2ax£ (3 + ecosx) 



~ . . „ 2Mx. . 

V v {x,P,e) = x + at (1 + ecosx), 

P 

f>/ % 2 c 2aMx 
V t {x,P,e) = at (l-ecosx) + 



£p 2 



P 



(1 + ecosx) 

J(x,P,e) = 1 - ^-(1 - ecosx) + ^r(l + ecosx) 2 
p p z 



2 ' 



(Al) 
(A2) 

(A3) 
(A4) 
(A5) 
(A6) 



The quantity x = x(a,p, e) [which also enters Eqs. ([3]) and ((U) for the particle's energy £ and angular- momentum C] 
is given by 



—N - sign(a)VN 2 ~ AFC 



2F 



1/2 



with 



(A7) 

>M(3 + e 2 )p 2 + A/ 2 (3 + e 2 ) 2 p-4Ma 2 (l-e 2 ) 2 ] , (A8) 

+ [Af 2 (3 + e 2 ) - a 2 ] p - Ma 2 (l + 3e 2 )} , (A9) 
) 2 . (A10) 

Appendix B: Comparison with previous results in the Schwarzschild case 

Canizares et al. (3(| recently presented numerical results from a calculation of the SSF using a pseudospectral 
collocation method. Their calculation was carried out in the time domain, restricting to the Schwarzschild case 
(a = 0). In Table HVl we show a comparison of their results with ours, showing a good agreement. 



F(p,e) 


1 






N(p,e) 


= 2 -\ 
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C(P) 


= (a* 
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Table IV: Comparison with Canizares et al. [3(| in the Schwarzs child case (a = 0). The SSF values are extracted at certain 
near-periastron points as specified in Table I of [3(J. Canizares et al. do not indicate error bars on their results; for our data 
all figures are significant. 
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